/*
 * The MIT License (MIT)
 *
 * Copyright (c) 2013 Jean-Sebastien Lerat (Jean-Sebastien.Lerat@ulb.ac.be)
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

/**
 *
 * @author Jean-Sebastien Lerat (Jean-Sebastien.Lerat@ulb.ac.be)
 * @version 1.00, 24/06/2013
 * Native implementation of the measure library for CoNet (http://psbweb05.psb.ugent.be/conet/)
 */

#include "measures.h"

MeasureFunction defineMeasure(JNIEnv * env,jobject config){
	jobject object	= getClass(env,config);

	jstring strName	= getStringName(env,object);
	const char *name= (*env)->GetStringUTFChars((JNIEnv *)env,strName,0);
	const char *simple_name = strrchr(name,'.')+1;
	
	MeasureFunction result;
	result.config	= NULL;
	
	if(       equals(simple_name,"BrayCurtis")){
		result.function =	&bray_curtis;
	} else if(equals(simple_name,"ChiSquare")){
		result.function =	&chi_square;
	} else if(equals(simple_name,"DistanceCorrelation")){
		result.function =	&distance_correlation;
	} else if(equals(simple_name,"EmpiricalEntropy")){
		const jfieldID field_ptr	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,config),"pseudoCount","D");
		const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptr);
		jdouble * ptr = (jdouble *)malloc(sizeof(jdouble));
		*ptr				= pseudoCount;
		result.config		= ptr;
		result.function =	&empirical_entropy;
	} else if(equals(simple_name,"Euclidean")){
		result.function =	&euclidean;
	} else if(equals(simple_name,"Hellinger")){
		result.function =	&hellinger;
	} else if(equals(simple_name,"HilbertSchmidt")){
		result.function =	&hilbert_schmidt;
	} else if(equals(simple_name,"JensenShannon")){
		const jfieldID field_ptr	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,config),"pseudoCount","D");
		const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptr);
		jdouble * ptr = (jdouble *)malloc(sizeof(jdouble));
		*ptr				= pseudoCount;
		result.config		= ptr;
		result.function =	&jensen_shannon;
	} else if(equals(simple_name,"Kendall")){
		result.function	= &kendall;
	} else if(equals(simple_name,"KullbackLeibler")){
		const jfieldID field_ptr	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,config),"pseudoCount","D");
		const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptr);
		jdouble * ptr = (jdouble *)malloc(sizeof(jdouble));
		*ptr				= pseudoCount;
		result.config		= ptr;
		result.function =	&kullback_leibler;
	} else if(equals(simple_name,"MIC")){
		result.function =	&mic;
	} else if(equals(simple_name,"MutualInformation")){
		const jfieldID field_ptr	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,config),"entropyEstimator","Lbe/ac/ulb/mlg/utils/measure/Entropy;");
		const jobject jEntropy	= (*env)->GetObjectField(env,config,(jfieldID)field_ptr);
		jobject jEntropyClass		= getClass(env,jEntropy);

		jstring strNameEntropy					= getStringName(env,jEntropyClass);
		const char *nameEntropy				= (*env)->GetStringUTFChars((JNIEnv *)env,strNameEntropy,0);
		const char *simple_nameEntropy	= strrchr(nameEntropy,'.')+1;

		void* ptr;
/* FIXME
		if(equals(simple_nameEntropy,"EmpiricalEntropy")){
			const jfieldID field_ptrPS	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,jEntropy),"pseudoCount","D");
			const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptrPS);
			ptr					= malloc(sizeof(Entropy)+sizeof(jdouble));
			*((*Entropy)ptr)	= &entropy_empirical;
			((*Entropy)ptr)++;
			*((jdouble*)ptr) = pseudoCount;
		}else if(equals(simple_nameEntropy,"SchurmannGrassbergerEntropy")){
			ptr = malloc(sizeof(Entropy));
			*((*Entropy)ptr) = &sg_entropy;
		}else if(equals(simple_nameEntropy,"ShannonEntropy")){
			const jfieldID field_ptrPS	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,jEntropy),"pseudoCount","D");
			const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptrPS);
			ptr					= malloc(sizeof(Entropy)+sizeof(jdouble));
			*((*Entropy)ptr)	= &entropy_shannon;
			((*Entropy)ptr)++;
			*((jdouble*)ptr) = pseudoCount;
		}
*/
		result.config		= ptr;
		result.function =	&mutual_information;
	} else if(equals(simple_name,"Pearson")){
		result.function =	&pearson;
	} else if(equals(simple_name,"PolynomialRegression")){
		result.function =	&polynomial_regression;
	} else if(equals(simple_name,"SchurmannGrassbergerEntropy")){
		result.function =	&sg_entropy;
	} else if(equals(simple_name,"ShannonEntropy")){
		const jfieldID field_ptr	= (*env)->GetFieldID(env,(*env)->GetObjectClass(env,config),"pseudoCount","D");
		const jdouble pseudoCount	= (*env)->GetDoubleField(env,config,(jfieldID)field_ptr);
		jdouble * ptr = (jdouble *)malloc(sizeof(jdouble));
		*ptr				= pseudoCount;
		result.config		= ptr;
//FIXME get pseudoCount from the EmpiricalEntropy class
		result.function =	&shannon_entropy;
	} else if(equals(simple_name,"Spearman")){
		result.function =	&spearman;
	}else if(equals(simple_name,"Steinhaus")){
		result.function =	&steinhaus;
	} else if(equals(simple_name,"VarianceOfLogRatios")){
		result.function =	&variance_log_ratios;
	} else {
		char exBuffer[256];
		(*env)->ReleaseStringUTFChars((JNIEnv *)env, strName, name);
		sprintf(exBuffer, "TODO appropriate error msg  %s %s print", "log", "or");
		(*env)->ThrowNew(env, (*env)->FindClass(env, "java/lang/Exception"), exBuffer);
		//return NULL;
	}
	(*env)->ReleaseStringUTFChars((JNIEnv *)env, strName, name);
	return result;
}

jdouble empirical_entropy(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize M){
	double frequencies[M];
	const jdouble* mat[2] = {A,B};
	int position[M];
	const long informations = countOccurenciesMat(mat, M, frequencies, M,position, computable);
	const int distinctValues = (int)(informations>>32);
	const int numberOfValues = (int)(informations>>32);

	return entropy_empirical(config,frequencies,M,distinctValues,numberOfValues);
}
jdouble shannon_entropy(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize M){
	double frequencies[M];
	const jdouble * mat[2] = {A,B};
	int position[M];
	const long informations = countOccurenciesMat(mat, M, frequencies, M,position, computable);
	const int distinctValues = (int)(informations>>32);
	const int numberOfValues = (int)(informations>>32);

	return entropy_shannon(config,frequencies,M,distinctValues,numberOfValues);
}
jdouble sg_entropy(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize M){
	double frequencies[M];
	const jdouble * mat[2] = {A,B};
	int position[M];
	const long informations = countOccurenciesMat(mat, M, frequencies, M,position, computable);
	const int distinctValues = (int)(informations>>32);
	const int numberOfValues = (int)(informations>>32);

	return entropy_sg(config,frequencies,M,distinctValues,numberOfValues);
}

jdouble polynomial_regression(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	//TODO not yet implemented
	return NAN;
}

jdouble distance_correlation(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	//TODO not yet implemented
	return NAN;
}

jdouble chi_square(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	//TODO not yet implemented
	return NAN;
}

jdouble mic(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	//TODO not yet implemented
	return NAN;
}

jdouble jensen_shannon(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble KL = 0.0,M;
	jsize index;
	jint count = N;
	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]){
				M	=   (A[index]+B[index])/2.0;
				KL	+= A[index]*log(A[index]/M + *((jdouble*)config))+B[index]*log(B[index]/M + *((jdouble*)config));
			}
	}else{
		for(index=0;index<N;index++){
			M	=   (A[index]+B[index])/2.0;
			KL	+= A[index]*log(A[index]/M + *((jdouble*)config))+B[index]*log(B[index]/M + *((jdouble*)config));
		}
	}
	//if(count < 1)return ERROR_VALUE;
	count<<=1; // uncomment for correct normalization
	//KL /= count;
	if(isinf(KL<0.0?-KL:KL)) return ERROR_VALUE;
	return KL;
}

jdouble hilbert_schmidt(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	int n = N;
	int i,j,k,l,m;

	double trace = 0.0;
	double tmp1,tmp2,tmp3;//used to "prevent" error during accumulations (large + small)

	double multA = 0.0 ,multB = 0.0;

	{
		double sumA = 0.0,sumB = 0.0;
		if(computable) {
			for(i=0;i<N;i++){
				if(!computable[i]){
					--n;
				}
				if(isnan(A[i])){
					sumA	+= A[i];
					multA	+= A[i]*A[i];
				}
				if(isnan(B[i])){
					sumB	+= B[i];
					multB	+= B[i]*B[i];
				}
			}
		}else{
			for(i=0;i<N;i++){
				sumA	+= A[i];
				multA	+= A[i]*A[i];
				sumB	+= B[i];
				multB	+= B[i]*B[i];
			}
		}
		multA = - 1.0  / (2.0 * ( multA - sumA*sumA )/n);
		multB = - 1.0  / (2.0 * ( multB - sumB*sumB )/n);
	}

	const double gamma0 = -1.0 / n;
	const double gamma1 = 1.0 + gamma0;

	if(computable) {
		int notNaN[n];
		for(i=0,j=0;i<n;i++,j++){
			while(!computable[j]) ++j;
			notNaN[i]=j;
		}
		for(i=0;i<n;i++){
			j = i;
			tmp3 = 0.0;
			for(m=0;m<n;m++){
				tmp2 = 0.0;
				for(l=0;l<n;l++){
					tmp1 = 0.0;
					for(k=0;k<n;k++){
						tmp1 += (k!=l?gamma0:gamma1) * GAUSSIAN_KERNEL(multA,A[notNaN[i]],A[notNaN[k]]);
					}
					tmp2 += tmp1 * GAUSSIAN_KERNEL(multB,B[notNaN[l]],B[notNaN[m]]);
				}
				tmp3 += tmp2 * (m!=i?gamma0:gamma1);
			}
			trace += tmp3;
		}
	}else{
		for(i=0;i<n;i++){
			j = i;
			tmp3 = 0.0;
			for(m=0;m<n;m++){
				tmp2 = 0.0;
				for(l=0;l<n;l++){
					tmp1 = 0.0;
					for(k=0;k<n;k++){
						tmp1 += (k!=l?gamma0:gamma1) * GAUSSIAN_KERNEL(multA,A[i],A[k]);
					}
					tmp2 += tmp1 * GAUSSIAN_KERNEL(multB,B[l],B[m]);
				}
				tmp3 += tmp2 * (m!=i?gamma0:gamma1);
			}
			trace += tmp3;
		}
	}
	return trace/(n*n);
}

jdouble variance_log_ratios(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	double mean		= 0.0;
	double square	= 0.0;
	double logRatio;
	int index;
	int count = N;

	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]){
				logRatio = log(A[index]/B[index]);//TODO slog ?
				mean	+= logRatio;
				square	+= logRatio*logRatio;
			}else count--;
	}else{
		for(index=0;index<N;index++){
			logRatio = log(A[index]/B[index]);
			mean	+= logRatio;
			square	+= logRatio*logRatio;
		}
	}

	if(count < 1)return ERROR_VALUE;

	mean/=count;
	square/=count;

	return 1.0 - exp(-ssqrt(square - mean*mean));
}

jdouble hellinger(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jsize index;

	jdouble hellinger = 0.0,diff;
	
	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]){
				diff		=  sqrt(A[index])-sqrt(B[index]);
				hellinger	+= diff*diff;
			}
	}else{
		for(index=0;index<N;index++){
			diff		=  sqrt(A[index])-sqrt(B[index]);
			hellinger	+= diff*diff;
		}
	}

	if(hellinger < 0.0) return ERROR_VALUE;

	return ssqrt(hellinger)*0.70710678118;//0.70710678118 = 1/sqrt(2)
}

jdouble euclidean(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble sum = 0.0,square;
	jsize i;
	
	if(computable != NULL){
		for(i=0;i<N;i++)
			if(computable[i]){
				square = A[i]-B[i];//faster than call power
				sum+=square*square;
			}
	}else{
		for(i=0;i<N;i++){
			square = A[i]-B[i];//faster than call power
			sum+=square*square;
		}
	}
	
	if(sum < 0.0) return ERROR_VALUE;

	return ssqrt(sum);
}

jdouble kendall(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble rankA[N];
	jdouble rankB[N];
	int rankJoint[N];

	const int n1 = naturalRank(computable,A,rankA,rankJoint,N);
	const int n2 = naturalRank(computable,B,rankB,rankJoint,N);

	int n0=0,n3=0,t3=0;

	const double upperBounds = N + 1;

	int index;

	if(computable != NULL){
		int previous;
		for(previous=0;previous < N && !computable[previous]; previous++){
			rankA[previous] = rankB[previous] = -INFINITY;
		}
		if(previous >= N) return ERROR_VALUE;//if previous is not < computable.length => all is NaN
		rankA[previous] = rankA[previous] * upperBounds  + rankB[previous];
		++n0;
		for(index=previous+1;index<N;index++)
			if(computable[index]){
				rankA[index] = rankA[index] * upperBounds  + rankB[index];
				if(abs(rankA[index]-rankA[previous])<1e-14){
					t3++;
				}else{
					n3 += t3*(t3+1);
					t3  = 0;
				}
				++n0;
				previous = index;
			}else{
				rankA[index] = rankB[index] = -INFINITY;
			}
	}else{
		rankA[0] = rankA[0] * upperBounds  + rankB[0];
		for(index=1;index<N;index++){
			rankA[index] = rankA[index] * upperBounds  + rankB[index];
			if(abs(rankA[index]-rankA[index-1])<1e-14){
				t3++;
			}else{
				t3  = 0;
			}
		}
		n0 = N;
	}
	n3 += t3*(t3+1);
	n3 >>= 1;
	n0 = (n0 * (n0-1)) >> 1;

	bubbleSort(rankA, rankB,N);

	const int swap = countSwap(rankB,N);

	//sort rankB (only indicies) based on the sorting of rankA (rankJoint)
	const double numerator		= n0 + ( -n1-n2+n3 - (swap<<1));//nc -nd = n0 - n1 - n2 + n3 - 2*S(y)
	const double denominator	= exp(0.5*(log(n0-n1)+log(n0-n2)));//overflow (Numerical Analysis course) (0.5 sqrt = exp(1/2))

	if(denominator<1e-16)return ERROR_VALUE;

	return numerator/denominator;
}

jdouble mutual_information(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	Entropy estimator		= config;
	const void* config2	= config++;

	long informations;
	const jsize joint_len = N;
	double frequencies[joint_len];
	int position[joint_len];
	const double * joint[2] = {(double*)A,(double*)B};

	informations = countOccurenciesVect(A, N,frequencies,joint_len,position,computable != NULL);
	const double entropyA = estimator(config2,frequencies,joint_len,(int)(informations>>32),(int)informations);

	informations = countOccurenciesVect(B, N,frequencies,joint_len,position,computable != NULL);
	const double entropyB = estimator(config2,frequencies,joint_len,(int)(informations>>32),(int)informations);

	informations = countOccurenciesMat(joint,N,frequencies,joint_len,position,computable);
	const double entropyAB = estimator(config2,frequencies,joint_len,(int)(informations>>32),(int)informations);

	return entropyA + entropyB - entropyAB;
}

jdouble bray_curtis(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble min = 0.0,sum = 0.0;
	jsize index;
	
	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]){
				min += fmin(A[index],B[index]);
				sum += A[index]+B[index];
			}
	}else{
		for(index=0;index<N;index++)
			min += fmin(A[index],B[index]);
			sum += A[index]+B[index];
	}

	if((sum>0?sum:-sum)<1e-16)return ERROR_VALUE;

	return 1.0 - 2.0 * min / sum;
}

jdouble steinhaus(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble min = 0.0,sum = 0.0;
	jsize index;

	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]){
				min += fmin(A[index],B[index]);
				sum += A[index]+B[index];
			}
	}else{
		for(index=0;index<N;index++)
			min += fmin(A[index],B[index]);
			sum += A[index]+B[index];
	}

	if((sum>0?sum:-sum)<1e-16)return ERROR_VALUE;

	return 2.0 * min / sum;
}

jdouble kullback_leibler(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble KL = 0.0;
	jsize index;
	jint count = N;
	if(computable != NULL){
		for(index=0;index<N;index++)
			if(computable[index]) KL += A[index]*log(A[index]/B[index]+ *((jdouble*)config))+B[index]*log(B[index]/A[index]+ *((jdouble*)config));
			else count--;
	}else{
		for(index=0;index<N;index++)
			KL += A[index]*log(A[index]/B[index]+ *((jdouble*)config))+B[index]*log(B[index]/A[index]+ *((jdouble*)config));
	}
	//if(count < 1)return ERROR_VALUE;
	count<<=1; // uncomment for correct normalization
	//KL /= count;
	if(isinf(KL<0.0?-KL:KL)) return ERROR_VALUE;
	return KL;
}

jdouble spearman(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	jdouble rankA[N];
	jdouble rankB[N];
	int positions[N];

	naturalRank(computable,A,rankA,positions,N);
	naturalRank(computable,B,rankB,positions,N);

	const jdouble avg_A	= fastMean(rankA,computable,N);
	const jdouble avg_B	= fastMean(rankB,computable,N);
	jsize index;
	
	//compute the squares
	jdouble square_B=0.0, mult_AB = 0.0, square_A = 0.0;
	jdouble tmp_A,tmp_B;


	if(computable != NULL){
		for(index=0; index<N; index++)
			if(computable[index]){
				tmp_A = rankA[index] - avg_A;
				tmp_B = rankB[index] - avg_B;

				square_A	+= tmp_A*tmp_A;
				square_B	+= tmp_B*tmp_B;
				mult_AB		+= tmp_A*tmp_B;
			}
	}else{
		for(index=0; index<N; index++){
			tmp_A = rankA[index] - avg_A;
			tmp_B = rankB[index] - avg_B;

			square_A	+= tmp_A*tmp_A;
			square_B	+= tmp_B*tmp_B;
			mult_AB		+= tmp_A*tmp_B;
		}
	}

	if((square_A>0?square_A:-square_A)<1e-16)return ERROR_VALUE;
	if((square_B>0?square_B:-square_B)<1e-16)return ERROR_VALUE;

	return mult_AB / ssqrt(square_A * square_B);
}

jdouble pearson(const void* config,const jdouble* A,const jdouble* B,const jboolean * computable,const jsize N){
	const jdouble avg_A	= fastMean(A,computable,N);
	const jdouble avg_B	= fastMean(B,computable,N);
	jsize index;
	
	//compute the squares
	jdouble square_B=0.0, mult_AB = 0.0, square_A = 0.0;
	jdouble tmp_A,tmp_B;
	
	if(computable != NULL){
		for(index=0; index<N; index++)
			if(computable[index]){
				tmp_A = A[index] - avg_A;
				tmp_B = B[index] - avg_B;
				
				square_A	+= tmp_A*tmp_A;
				square_B	+= tmp_B*tmp_B;
				mult_AB	+= tmp_A * tmp_B;
			}
	}else{
		for(index=0; index<N; index++){
			tmp_A = A[index] - avg_A;
			tmp_B = B[index] - avg_B;
			
			square_A	+= tmp_A*tmp_A;
			square_B	+= tmp_B*tmp_B;
			mult_AB	+= tmp_A * tmp_B;
		}
	}
	return mult_AB / exp(0.5*log(square_A*square_B));
}

void bubbleSortPositional(jdouble * comparable,int * toSort,const jsize len){
	double tmpDouble;
	int		end		= len,
			index,
			currentEnd,
			tmpInt;

	do{
		currentEnd = 0;
		for(index = 1;index<end;index++)
			if(comparable[index-1] > comparable[index]){
				tmpDouble			= comparable[index-1];
				comparable[index-1]	= comparable[index];
				comparable[index]	= tmpDouble;

				tmpInt			= toSort[index-1];
				toSort[index-1]	= toSort[index];
				toSort[index]	= tmpInt;

				currentEnd = index;
			}
		end = currentEnd;
	}while(end>0);
}

void bubbleSort(jdouble * comparable,jdouble * toSort,const jsize len){
	double tmpDouble;
	int		end		= len,
			index,
			currentEnd,
			tmpInt;

	do{
		currentEnd = 0;
		for(index = 1;index<end;index++)
			if(comparable[index-1] > comparable[index]){
				tmpDouble			= comparable[index-1];
				comparable[index-1]	= comparable[index];
				comparable[index]	= tmpDouble;

				tmpDouble		= toSort[index-1];
				toSort[index-1]	= toSort[index];
				toSort[index]		= tmpDouble;

				currentEnd = index;
			}
		end = currentEnd;
	}while(end>0);
}

int countSwap(jdouble * comparable,const jsize len){
	double tmpDouble;
	int		end		= len,
			index,
			currentEnd,
			swaps = 0;

	do{
		currentEnd = 0;
		for(index = 1;index<end;index++)
			if(comparable[index-1] > comparable[index]){
				swaps++;
				tmpDouble			= comparable[index-1];
				comparable[index-1]	= comparable[index];
				comparable[index]	= tmpDouble;

				currentEnd = index;
			}
		end = currentEnd;
	}while(end>0);
	return swaps;
}

int mergesortCounting(const jboolean * computable,int * toSort,jdouble * comparable,const jsize len){
	int buff[len];
	int * swapper = buff;

	int tapeSize,swapCount = 0;
	int index,leftTape,rightTape,end;
	int ptrLeftTape,ptrRightTape,position;

	for(tapeSize = 1; tapeSize < len ; tapeSize <<= 1){
		for(index = 0; index < len ; index += tapeSize<<1){
			leftTape	= index;
			rightTape	= index+tapeSize;
			end			= rightTape+tapeSize;
			if(rightTape > len) rightTape = len;
			if(end > len) end = len;

			ptrLeftTape		= leftTape;
			ptrRightTape	= rightTape;

			for(position = leftTape; position < end; position++){
				if(ptrLeftTape >= rightTape || (ptrRightTape < end && comparable[toSort[ptrLeftTape]] > comparable[toSort[ptrRightTape]])){
					swapper[position]	=  toSort[ptrRightTape++];
					swapCount			+= rightTape - ptrLeftTape;
				}else{
					swapper[position] = toSort[ptrLeftTape++];
				}
			}

		}

		int * tmp	= toSort;
		toSort		= swapper;
		swapper	= toSort;
	}

	return swapCount;
}

int naturalRank(const jboolean * computable,const jdouble * vector,jdouble * ranks,int * position,const jsize len){
	int index,count=0,ties=0;

	if(computable != NULL){
		for(index = 0;index < len; index++) {
			position[index]	= index;
			ranks[index]	= computable[index]?vector[index]:INFINITY;
		}
	}else{
		for(index = 0;index < len; index++) {
			position[index]	= index;
			ranks[index]	= vector[index];
		}
	}

	bubbleSortPositional(ranks,position,len);

	int pos = 1, previous = 0;

	if(computable != NULL){
		int startTie;
		while(!computable[position[previous]]) previous++;
		startTie = previous;

		ranks[position[previous]] = pos;

		for(index=previous+1; index<len; index++){
			if(computable[position[index]]){
				if(vector[position[index]] > vector[position[previous]]){// different but > it's ok because it's sorted
					pos++;
					if(count){
						pos+= count;
						const jdouble average =  (2.0 * ranks[position[previous]] + count) * 0.5;
						previous = 0;
						count++;
						for(count=startTie;count<index ; count++){
							if(computable[position[count]]){
								ranks[position[count]] = average;
								previous++;
							}
						}
						ties	+= (previous*(previous-1));
						count	=  0;
					}
				} else {
					count++;
				}
				ranks[position[index]]	= pos;
				previous					= index;
			}
		}
		if(count){
			pos+= count;
			const jdouble average =  (2.0 * ranks[position[previous]] + count) * 0.5;

			previous = 0;
			count++;
			for(count=startTie;count<index ; count++){
				if(computable[position[count]]){
					ranks[position[count]] = average;
					previous++;
				}
			}
			ties	+= (previous*(previous-1));
			count	=  0;
		}
	}else{
		ranks[position[0]] = pos;
		for(index=1; index<len; index++) {
			if(vector[position[index]] > vector[position[index-1]]){// different but > it's ok because it's sorted
				pos = index + 1;
				if(count){
					const jdouble eq_value =  (2.0 * ranks[position[index-1]] + count) * 0.5;
					count++;
					ties += count*(count-1);
					while(count){
						count--;
						ranks[position[index-count]] = eq_value;
					}
				}
			} else {
				count++;
			}
			ranks[position[index]] = pos;
		}

		if(count){
			const jdouble eq_value =  (2.0 * ranks[position[index-1]] + count) * 0.5;
			count++;
			ties += count*(count-1);
			while(count){
				count--;
				ranks[position[index-count]] = eq_value;
			}
		}
	}
	return ties>>1;
}


jdouble fastSum(const jdouble* vector,const jboolean * computable,jsize len){
	jdouble sum = 0.0;
	if(computable != NULL){
		vector += len-1;
		computable += len-1;
		while(len--){
			if(*computable) sum += (*vector);
			vector--;
			computable--;
		}
	}else{
		while(len--){
			sum += (*vector);
			vector++;
		}
	}
	return sum;
}

jdouble fastMean(const jdouble* vector,const jboolean * computable,jsize len){
	jdouble sum = 0.0;
	jint count = len;
	if(computable != NULL){
		vector += len-1;
		computable += len-1;
		while(len--){
			if(*computable)sum += (*vector);
			else count--;
			vector--;
			computable--;
		}
	}else{
		while(len--){
			sum += (*vector);
			vector++;
		}
	}
	return sum/count;
}

///////////////////////////////////
long countOccurenciesMat(const jdouble** vector,const jsize M,double * frequencies,const jsize freq_len,int * position,const jboolean * computable){//TODO update with java version
	jdouble jointValues[M];
	{
		int j;
		if(computable){
			for(j=0;j<M;j++){
				jointValues[j] = computable[j] ? (jdouble) ((((long)vector[0][j])<<32) | ((long)vector[1][j])) : (jdouble) NAN;
			}
		}else{
			for(j=0;j<M;j++)
				jointValues[j] = (jdouble) ((((long)vector[0][j])<<32) | ((long)vector[1][j]));
		}
	}
	return countOccurenciesVect((const jdouble*)jointValues, M,frequencies,M, position, computable != NULL);
}

long countOccurenciesVect(const jdouble* vector,const jsize len,double * frequencies,const jsize freq_len,int * position,const jboolean hasMissing){
	int index,distinct = 0,size = len;

	if(hasMissing){
		for(index = 0;index < len; index++) {
			position[index]	= index;
			if(isnan(vector[index])){
				frequencies[index] = INFINITY;
				--size;
			}else{
				frequencies[index]	= vector[index];
			}
		}
	}else{
		for(index = 0;index < len; index++) {
			position[index]		= index;
			frequencies[index]	= vector[index];
		}
	}

	bubbleSortPositional(frequencies,position,len);
	for(index = 0;index < freq_len; index++){
		frequencies[index] = 0.0;
	}

	int count = 1;
	int sum = 0;
	for(index=1;index<size;index++){
		if(abs(vector[position[index-1]]-vector[position[index]])<EPSILON){//it is the same
			++count;
		}else{
			++distinct;
			sum += count;
			frequencies[position[index-count]] = count;
			count = 1;
		}
	}
	++distinct;
	sum += count;
	frequencies[position[index-count]] = count;

	return (((long)distinct) << 32) | sum;
}

/*
 * digamma function inspired from commons-apache
 */
double digamma(double x) {
	if (x <= 1E-5)
		return -EULER_MASCHERONI - 1.0 / x;

	double acc = 0.0;
	while(x < 49){
		acc += 1.0 / x;
		x	+= 1.0;
	}

	const double inv = 1.0 / (x * x);
	return log(x) - 0.5 / x - inv * ((1.0 / 12.0) + inv * (1.0 / 120.0 - inv / 252.0)) - acc;
}

double entropy_sg(const void* config,double *frequencies,const jsize len,const int distinctValues,const int numberOfValues){
	const double pseudoCount			= 1.0 / distinctValues;
	const double binsXpseudoCount	= distinctValues * pseudoCount;

	double sum = 0.0;
	int i;
	for(i=0;i<len;i++){
		if(frequencies[i]>0.0)
			sum += (frequencies[i] + pseudoCount) * (digamma( numberOfValues + binsXpseudoCount + 1.0) - digamma(frequencies[i] + pseudoCount + 1.0));
	}

	return sum / (numberOfValues + binsXpseudoCount);
}
double entropy_shannon(const void* config,double *frequencies,const jsize len,const int distinctValues,const int numberOfValues){
	return entropy_shannon(config,frequencies,len,distinctValues,numberOfValues)/LOG_E_2;
}
double entropy_empirical(const void* config,double *frequencies,const jsize len,const int distinctValues,const int numberOfValues){
	double entropy = 0.0;
	int i;
	double freq;
	for(i=0;i<len;i++){
		freq	= (frequencies[i]+ *((jdouble*)config))/numberOfValues;
		if(freq > 0.0)
			entropy += freq * log(freq);
	}
	return -entropy;
}
